additional_analyses

Author

Liz Loutrage

Data preparation

Data summary

Traits correlation

  • remove the traits that are too correlated?
Code
M <-cor(numeric_traits[, c(-1)])

ggcorrplot::ggcorrplot(M, hc.order = TRUE, type = "lower",
                       lab = TRUE, tl.cex = 9, lab_size = 3)

  • Traits types

CWM bootstrap

By depth layer

  • each trawl is a replica

  • Utilisation traitstrap R package

  • Bootstrap non parametric, échantillonnage dans les proportions de la biomasse des espèces

  • pour le remplissage des valeurs manquantes : si manque données pour 1 individus alors va chercher la valeurs sur un ind de la même profondeur, si pas possible alors ind à une autre profondeur mais de la même couche et si impossible ind d’une autre couche de profondeur

Code
spxcom.matrix <-  depth_fish_biomass %>% 
  t() %>% 
  as.data.frame() %>% 
  relocate("Epipelagic", "Upper mesopelagic", "Lower mesopelagic","Bathypelagic" ) %>% 
  tibble::rownames_to_column("species") %>% 
  arrange(species) %>% 
  tibble::column_to_rownames("species") %>% 
  as.matrix()

spxtraits.matrix <- fish_traits %>%
  mutate(across(16:23, ~ case_when(. == "P" ~ 1, 
                                   . == "A" ~ 0, 
                                   TRUE ~ as.numeric(.))),
         across(24, ~ case_when(. == "A" ~ 1, 
                                . == "B" ~ 2, 
                                . == "C" ~ 3, 
                                TRUE ~ as.numeric(.)))) %>% 
  select(-c(gill_raker_types, oral_gape_axis)) %>% 
  as.matrix()

# Remove the "[,1]" suffix from column names
names(spxtraits.matrix) <- gsub("[,1]", "", names(spxtraits.matrix))

#check rownames
#rownames(spxtraits.matrix) == rownames(spxcom.matrix)

result_CWM <- FD::functcomp(spxtraits.matrix, t(spxcom.matrix)) 
#FD::functcomp(spxtraits.matrix, t(spxcom.matrix), CWM.type = "all")

#  Calculate Total biomass
total_biomass <- colSums(spxcom.matrix)

#  Calculate Relative biomass
sp_rel_biomass <- t(spxcom.matrix) / total_biomass

# Transpose the Relative biomass Matrix for Display
t_sp_rel_biomass <- t(sp_rel_biomass)

total_sum <- colSums(t(sp_rel_biomass))

# Initialize an empty data frame to store results
CWM_df <- data.frame(
  depth_layer = character(),
  trait = character(),
  total_sum = numeric(),
  weighted_mean = numeric(),
  stringsAsFactors = FALSE
)

library(traitstrap)

# Trait 
trait_boot <- morpho_data%>% 
  inner_join(metadata) %>% 
  select(-c(individual_code, years, longitude_start,
            latitude_start, longitude_end, longitude_end,
            volume_filtered, distance)) %>% 
  mutate(
    eye_size = eye_diameter / head_depth,
    orbital_length = eye_diameter / standard_length,
    oral_gape_surface = mouth_width * mouth_depth / body_width * body_depth,
    oral_gape_shape = mouth_depth / mouth_width,
    oral_gape_position = distance_upper_jaws_bottom_head / head_depth,
    lower_jaw_length = lower_jaw_length / standard_length,
    head_length = head_length / standard_length,
    body_depth = body_depth / standard_length,
    pectoral_fin_position = distance_pectoral_bottom_body / body_depth_pectoral_insertion,
    pectoral_fin_insertion = prepectoral_length / standard_length,
    transversal_shape = body_depth / body_width,
    dorsal_fin_insertion = predorsal_length / standard_length,
    eye_position = eye_height / head_depth,
    operculum_volume = operculum_depth / operculum_width,
    gill_outflow = operculum_width,
    caudal_throttle_width = caudal_peduncle_min_depth
  ) %>%
  select(
    depth,
    species,
    eye_size,
    orbital_length,
    gill_outflow,
    oral_gape_surface,
    oral_gape_shape,
    oral_gape_position,
    lower_jaw_length,
    head_length,
    body_depth,
    pectoral_fin_position,
    pectoral_fin_insertion,
    transversal_shape,
    caudal_throttle_width,
    dorsal_fin_insertion,
    eye_position,
    operculum_volume,
    ventral_photophores, 
    gland_head,
    chin_barbel, 
    small_teeth, 
    large_teeth, 
    fang_teeth, 
    retractable_teeth, 
    internal_teeth
  ) %>%
  mutate_at(vars(ventral_photophores, 
                 gland_head,
                 chin_barbel, 
                 small_teeth, 
                 large_teeth, 
                 fang_teeth, 
                 retractable_teeth, 
                 internal_teeth), 
            funs(ifelse(. == "P", 1, ifelse(. == "A", 0, .)))) %>% 
  mutate(across(all_of(c("ventral_photophores", 
                         "gland_head",
                         "chin_barbel", 
                         "small_teeth", 
                         "large_teeth", 
                         "fang_teeth", 
                         "retractable_teeth", 
                         "internal_teeth")), as.numeric)) %>% 
  tidyr::pivot_longer(!c(species,depth), names_to = "trait", values_to = "values")%>%
  # add column with depth layer
  mutate(
    depth_layer = case_when(
      between(depth, 0, 174) ~ "Epipelagic",
      between(depth, 175, 699) ~ "Upper mesopelagic",
      between(depth, 700, 999) ~ "Lower mesopelagic",
      between(depth, 1000, 2000) ~ "Bathypelagic"))

# Community 
community <-  rbind(data_biomass_2002_2019, data_biomass_2021, data_biomass_2022)%>%
  as.data.frame()%>%
  rename("species"="Nom_Scientifique",
         "station"="Code_Station") %>%
  left_join(metadata) %>%
  select(species, Tot_V_HV, depth, volume_filtered)%>%
  # add column with depth layer
  mutate(
    depth_layer = case_when(
      between(depth, 0, 174) ~ "Epipelagic",
      between(depth, 175, 699) ~ "Upper mesopelagic",
      between(depth, 700, 999) ~ "Lower mesopelagic",
      between(depth, 1000, 2000) ~ "Bathypelagic"))%>%
  replace(is.na(.), 0)%>%
  group_by(species, depth)%>%
  mutate(biomass=sum(Tot_V_HV))%>%
  select(-c(Tot_V_HV))%>%
  distinct()%>%
  select(-c(volume_filtered)) %>% 
  filter(biomass>0) %>% 
  mutate(species = case_when(
    species == "Nannobrachium_atrum"~"Lampanyctus_ater",
    species == "Cyclothone"~"Cyclothone_sp",
    species == "Stomias_boa_boa"~"Stomias_boa",
    TRUE ~ species
  )) 

trait_filling <- trait_fill(
  # input data (mandatory)
  comm = community,
  traits = trait_boot,
  
  # specifies columns in your data (mandatory)
  abundance_col = "biomass",
  taxon_col = "species",
  trait_col = "trait",
  value_col = "values",
  
  # specifies sampling hierarchy
  scale_hierarchy = c("depth_layer", "depth"),
  
  # min number of samples
  min_n_in_sample = 5
)

# run nonparametric bootstrapping
np_bootstrapped_moments <- trait_np_bootstrap(
  trait_filling, 
  nrep = 100
)

np_bootstrapped_moments_plot <- np_bootstrapped_moments %>% 
  mutate(trait= gsub("_"," ", trait))%>%
  filter(
    trait %in% c(
      "caudal throttle width",
      "oral gape surface",
       "gill outflow",
      "large teeth",
       "small teeth",
      "orbital length",  
      "eye size",
      "transversal shape",
      "operculum volume",
      "ventral photophores",
      "internal teeth",
      "eye position",
      "oral gape shape"
    )
  )

np_bootstrapped_moments_plot$trait <- factor(
  np_bootstrapped_moments_plot$trait,
  levels = c(
    "caudal throttle width",
     "oral gape surface",
    "gill outflow",
    "large teeth",
    "small teeth",
    "orbital length",
    "eye size",
    "transversal shape",
    "operculum volume",
    "ventral photophores",
    "internal teeth",
    "eye position",
    "oral gape shape"
  )
)

np_bootstrapped_moments_plot$depth_layer <- factor(np_bootstrapped_moments_plot$depth_layer, 
                                                   levels = c("Epipelagic", "Upper mesopelagic", "Lower mesopelagic", "Bathypelagic"))

# Mean
ggplot(np_bootstrapped_moments_plot, aes(x=depth_layer, y=mean)) +
  geom_boxplot(aes(col=depth_layer, fill=depth_layer), alpha=0.08) +
  scale_fill_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
  scale_color_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
  guides(col="none", fill="none")+
  labs(y="Bootstrapped CWM")+
  facet_wrap(~trait, scales = "free")+
  theme_light()+
  theme(axis.text.x = element_blank(), 
        axis.title.x = element_blank(), 
        panel.grid.minor =  element_blank(),
        axis.title.y.left =  element_text(size =14), 
        strip.text = element_text( face="bold"),  
        axis.title.y = element_text(size=15),
        axis.text.y = element_text(size=14),
        strip.background=element_rect(fill="white"),
        strip.text.x = element_text(size = 14, color = "black"))

Code
ggsave("CWM_boot.png", path = "figures/additional_analyses", dpi = 700, height = 10, width = 12)

PCA CWM

Code
sum_boot_moment <- trait_summarise_boot_moments(
  np_bootstrapped_moments
) %>% 
  ungroup() %>% 
  filter(
    trait %in% c(
      "caudal_throttle width",
      "oral_gape_surface",
      "gill_outflow",
      "large_teeth",
      "small_teeth",
      "orbital_length",  
      "eye_size",
      "transversal_shape",
      "operculum_volume",
      "ventral_photophores",
      "internal_teeth",
      "eye_position",
      "oral_gape_shape"
    )
  ) %>% 
  mutate(trait= gsub("_"," ", trait)) %>% 
  select(c(trait, depth_layer, mean)) %>% 
  group_by(
    depth_layer, trait
  ) %>% 
  summarise(mean_value= mean(mean)) %>% 
  distinct() %>% 
  ungroup() %>% 
  tidyr::pivot_wider(id_cols=depth_layer, values_from = "mean_value", names_from = "trait") %>% 
  tibble::column_to_rownames(var = "depth_layer")

res.pca <- FactoMineR::PCA(sum_boot_moment, graph = FALSE)

factoextra::fviz_pca_biplot(res.pca, repel = TRUE,
                            col.var = "gray50", 
                            col.ind = "black", 
                            pointsize = 2,
                            arrowsize = 0.2,
                            title = "", 
                            label = "var")+
  theme_light()+
  theme(panel.grid.minor =  element_blank(),
        panel.grid.major = element_blank())

Code
ggsave("PCA_CWM.png", path = "figures/additional_analyses", dpi = 700, height = 5, width = 7)

Summarizes bootstrapping output

Extracting raw distributions

Code
# run nonparametric bootstrapping
raw_dist_np <- trait_np_bootstrap(
  filled_traits = trait_filling,
  raw = TRUE
)

raw_dist_np$depth_layer <- factor(
  raw_dist_np$depth_layer,
  levels = c(
    "Epipelagic",
    "Upper mesopelagic",
    "Lower mesopelagic",
    "Bathypelagic"
  )
) 

raw_dist_np <- raw_dist_np%>% 
    mutate(trait= gsub("_"," ", trait)) %>% 
  filter(!trait%in% c("chin barbel", "fang teeth", 
                      "gland head", "internal teeth",
                      "large teeth", "retractable teeth", 
                      "small teeth", "ventral photophores"))

ggplot(raw_dist_np, aes(x = log(values), fill = depth_layer, col= depth_layer)) +
  geom_density(alpha = 0.1) +
  scale_color_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
  scale_fill_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
  labs(x = " log(Trait value)") +
  guides(col="none", fill="none")+
  facet_wrap(facets = vars(trait), scales = "free")+
  theme_light()+
  theme(strip.text.x = element_text(size = 10, face = "bold", color="black"),
    strip.background = element_rect(fill = "white"),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 10))

Functional redundancy

  • Species Richness (N)

  • Simpson diversity (D) : probability that two individuals randomly selected from a sample will belong to the same species

  • Rao diversity (Q) : Sum of distances between pairs of randomly chosen species in trait space weighted by a relative abundance

  • Functionnal redundancy Ricotta et al. (2016) = 1-Q/D

Code
#Depth biomass data 
depth_fish_biomass_indices <- rbind(data_biomass_2002_2019, data_biomass_2021, data_biomass_2022)%>%
  as.data.frame()%>%
  rename("species"="Nom_Scientifique",
         "station"="Code_Station") %>%
  left_join(metadata) %>%
  select(species, Tot_V_HV, depth, volume_filtered, depth)%>%
  # divise biomass by the volume filtered at each trawl (g.m3)
  mutate(biomass_cpu=(Tot_V_HV/volume_filtered)*1000)%>%
  select(species, depth, biomass_cpu)%>%
  replace(is.na(.), 0)%>%
  group_by(species, depth)%>%
  mutate(biomass=sum(biomass_cpu))%>%
  select(-c(biomass_cpu))%>%
  distinct()%>%
  arrange(depth) %>% 
  tidyr::pivot_wider(names_from = species, values_from = biomass)%>%
  replace(is.na(.), 0)%>%
  tibble::column_to_rownames(var = "depth")%>%
  #change species name
  rename("Lampanyctus_ater"="Nannobrachium_atrum")%>%
  rename("Cyclothone_sp"="Cyclothone")%>%
  rename("Stomias_boa"="Stomias_boa_boa") %>%
  as.matrix()

# Species richness ----
nsp <- as.data.frame(vegan::specnumber(depth_fish_biomass_indices)) %>%
  tibble::rownames_to_column(var = "depth") %>% 
  mutate(depth=as.numeric(depth))
colnames(nsp) <- c("depth", "species_richness")

# Diversity indices
diversity_indices_list <- SYNCSA::rao.diversity(depth_fish_biomass_indices, traits = fish_traits) 

diversity_indices <- data.frame(diversity_indices_list[c(2:3)]) %>% 
  tibble::rownames_to_column(var = "depth") %>% 
  mutate(depth=as.numeric(depth)) %>% 
  inner_join(nsp) %>% 
  mutate(functionnal_redundancy_ricotta= 1-(FunRao/Simpson)) %>% 
  tidyr::pivot_longer(! depth, names_to = "indices") %>% 
  # add column with depth layer
  mutate(
    depth_layer = case_when(
      between(depth, 0, 174) ~ "Epipelagic",
      between(depth, 175, 699) ~ "Upper mesopelagic",
      between(depth, 700, 999) ~ "Lower mesopelagic",
      between(depth, 1000, 2000) ~ "Bathypelagic"))

diversity_indices$depth_layer <-
  factor(
    diversity_indices$depth_layer,
    levels = c(
      "Epipelagic",
      "Upper mesopelagic",
      "Lower mesopelagic",
      "Bathypelagic"
    )
  )

diversity_indices$indices <- factor(
  diversity_indices$indices,
  levels = c(
    "species_richness",
    "Simpson",
    "FunRao",
    "functionnal_redundancy_ricotta"
  ),
  labels = c(
    "Species richness (N)",
    "Simpson diversity (D)",
    "Rao diversity (Q)",
    "Functional Redundancy R=1-(Q/D)"
  )
)

ggplot(diversity_indices, aes(x = depth, y = value)) +
  geom_point(alpha = 0.5, size = 1.5) +
  geom_smooth(method = "lm", col = "#008E9B", alpha = 0.1) +
  facet_wrap(~indices, scales = "free") +
  ggpmisc::stat_poly_eq(formula = y ~ x, 
                        aes(label = paste(..rr.label.., ..p.value.label.., ..n.label.., sep = "*`,`~")), 
                        parse = TRUE, 
                        hjust = 0, vjust = 1) +
  theme_light() +
  theme(axis.title.y.left = element_text(size = 14), 
        strip.text = element_text(size = 14, face = "bold"),  
        legend.title = element_text(size = 11),  
        legend.text = element_text(size = 11), 
        axis.title.y = element_text(size = 11),
        axis.text.y = element_text(size = 12),
        strip.background = element_rect(fill = "white"),
        strip.text.x = element_text(size = 11, face = "bold", color = "black"))

Code
ggplot(diversity_indices, aes(x =depth_layer, y=value))+
  geom_boxplot(aes(col=depth_layer, fill=depth_layer), alpha=0.1) +
  scale_fill_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
  scale_color_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A"))+
  facet_wrap(~indices, scales = "free")+
  theme_light()+
  guides(col="none", fill="none")+
  theme(axis.text.x = element_blank(), 
        axis.title.x = element_blank(), 
        axis.title.y.left =  element_text(size =14), 
        strip.text = element_text(size =14, face="bold"),  
        legend.title = element_text(size =11),  
        legend.text = element_text(size =11), 
        axis.title.y = element_text(size=11),
        axis.text.y = element_text(size=12),
        strip.background=element_rect(fill="white"),
        strip.text.x = element_text(size = 10, face = "bold", color = "black"))

Code
# # Data Preparation 
# station_morpho <- morpho_data %>% 
#   select(station) %>% 
#   filter(!station%in%c("F0315", "F0367","X0568")) %>%  
#   arrange(station) %>%
#   distinct()
# 
# station_sp <- rbind(data_biomass_2002_2019, data_biomass_2021, data_biomass_2022) %>%
#   as.data.frame() %>%
#   rename("species"="Nom_Scientifique",
#          "station"="Code_Station") %>%
#   left_join(metadata) %>%
#   select(species, Tot_V_HV, volume_filtered, station) %>%
#   # divise biomass by the volume filtered at each trawl (g.m3)
#   mutate(biomass_cpu=(Tot_V_HV/volume_filtered)*1000) %>%
#   select(species, biomass_cpu, station) %>%
#   replace(is.na(.), 0) %>%
#   group_by(species, station) %>%
#   mutate(biomass=sum(biomass_cpu)) %>%
#   select(-c(biomass_cpu)) %>%
#   distinct() %>%
#   tidyr::pivot_wider(names_from = species, values_from = biomass) %>%
#   replace(is.na(.), 0) %>%
#   arrange(station) %>%
#   filter(station %in% station_morpho$station) %>%
#   tibble::column_to_rownames(var = "station") %>%
#   #change species name
#   rename("Lampanyctus_ater"="Nannobrachium_atrum") %>%
#   rename("Cyclothone_sp"="Cyclothone") %>%
#   rename("Stomias_boa"="Stomias_boa_boa") %>%
#   as.matrix()
# 
# # Null Model Analysis
# nb_rep <- 1000
# 
# resultsRandom_FunRao <-
#   matrix(
#     NA,
#     nrow = nrow(station_sp),
#     ncol = nb_rep,
#     dimnames = list(rownames(station_sp),
#                     paste0("Sim.", 1:nb_rep))
#   )
# 
# resultsRandom_FunRedundancy <-
#   matrix(
#     NA,
#     nrow = nrow(station_sp),
#     ncol = nb_rep,
#     dimnames = list(rownames(station_sp),
#                     paste0("Sim.", 1:nb_rep))
#   )
# 
# for (rep in 1:nb_rep) {
#   randomize_mx <- picante::randomizeMatrix(samp = station_sp,
#                                            null.model =  "frequency",
#                                            iterations = 100)
#   
#   sim_cal <- SYNCSA::rao.diversity(randomize_mx , traits = fish_traits)
#   sim_FunRao <- sim_cal$FunRao
#   sim_FunRedundancy <- sim_cal$FunRedundancy
#   
#   resultsRandom_FunRao[, rep] <- sim_FunRao
#   resultsRandom_FunRedundancy[, rep] <- sim_FunRedundancy
# }
# 
# # Comparison with Observed Values
# obs <- SYNCSA::rao.diversity(station_sp, fish_traits)
# obs_FunRao <- obs$FunRao
# obs_FunRedundancy <- obs$FunRedundancy
# 
# # Calculate mean and standard deviation of null model FD values for FunRao
# meanNull_FunRao <- rowMeans(resultsRandom_FunRao, na.rm = TRUE)
# sdNull_FunRao <- apply(resultsRandom_FunRao, 1, sd, na.rm = TRUE)
# 
# # Calculate mean and standard deviation of null model FD values for FunRedundancy
# meanNull_FunRedundancy <- rowMeans(resultsRandom_FunRedundancy, na.rm = TRUE)
# sdNull_FunRedundancy <- apply(resultsRandom_FunRedundancy, 1, sd, na.rm = TRUE)
# 
# # Calculate effect size and standardized effect size for FunRao
# ES_FunRao <- obs_FunRao - meanNull_FunRao
# SES_FunRao <- ES_FunRao / sdNull_FunRao
# 
# # Calculate effect size and standardized effect size for FunRedundancy
# ES_FunRedundancy <- obs_FunRedundancy - meanNull_FunRedundancy
# SES_FunRedundancy <- ES_FunRedundancy / sdNull_FunRedundancy
# 
# # Combine the results into a data frame
# result_fd_df <- data.frame(
#   station = rownames(station_sp),
#   SES_FunRao = SES_FunRao,
#   SES_FunRedundancy = SES_FunRedundancy
# ) %>%
#   inner_join(metadata) %>%
#   select(
#     -c(
#       latitude_start,
#       longitude_start,
#       latitude_end,
#       longitude_end,
#       duration,
#       distance,
#       volume_filtered,
#       years,
#       comments,
#       station
#     )
#   ) %>% 
#   mutate(
#     depth_layer = case_when(
#       between(depth, 0, 174) ~ "Epipelagic",
#       between(depth, 175, 699) ~ "Upper mesopelagic",
#       between(depth, 700, 999) ~ "Lower mesopelagic",
#       between(depth, 1000, 2000) ~ "Bathypelagic"
#     )
#   ) %>% 
#   tidyr::pivot_longer(!c(depth, depth_layer), values_to = "values", names_to = "indice")
# 
# result_fd_df$depth_layer <- factor(result_fd_df$depth_layer, 
#                                    levels = c("Epipelagic", "Upper mesopelagic",
#                                               "Lower mesopelagic", "Bathypelagic"))
# 
# # Plot
# ggplot(result_fd_df, aes(x = depth_layer)) +
#   geom_point(aes(y = values, col=depth_layer), alpha=0.4) +
#   facet_wrap(~indice)+
#   geom_boxplot(aes(y = values, col=depth_layer, fill=depth_layer), alpha=0.07) +
#   scale_color_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A")) +
#   scale_fill_manual(values = c("#FEA520","#D62246","#6255B4","#3C685A")) +
#   labs(x="") +
#    theme_light()+
#   guides(col="none", fill="none")+
#   theme(axis.text.x = element_blank(), 
#         axis.title.x = element_blank(), 
#         axis.title.y.left =  element_text(size =14), 
#         strip.text = element_text(size =14, face="bold"),  
#         legend.title = element_text(size =11),  
#         legend.text = element_text(size =11), 
#         axis.title.y = element_text(size=11),
#         axis.text.y = element_text(size=12),
#         strip.background=element_rect(fill="white"),
#         strip.text.x = element_text(size = 10, face = "bold", color = "black"))+
#   geom_hline(yintercept=0, linetype="dashed", color = "black", linewidth=0.8)

SES Functional diversity indices

  • Standard Effect Size (SES): to eliminate the influence of species richness on the functional diversity indices (Mouchet et al., 2010). Measures the deviation from the random expectation in standard deviation units

  • null model frequency: Randomize community data matrix abundances (here biomasss) within species (maintains species occurence frequency)

  • FRic Functional Richness: the proportion of functional space filled by species of the studied assemblage, i.e. the volume inside the convex-hull shaping species. To compute FRic the number of species must be at least higher than the number of functional axis + 1.

  • FDis Functional Dispersion: the biomass weighted deviation of species traits values from the center of the functional space filled by the assemblage i.e. the biomass-weighted mean distance to the biomass-weighted mean trait values of the assemblage.

  • FDiv Functional Divergence: the proportion of the biomass supported by the species with the most extreme functional traits i.e. the ones located close to the edge of the convex-hull filled by the assemblage.

  • FEve Functional Evenness: the regularity of biomass distribution in the functional space using the Minimum Spanning Tree linking all species present in the assemblage.

  • FOri Functional Originality: the weighted mean distance to the nearest species from the global species pool.

  • FSpe Functional Specialization: the biomass weighted mean distance to the mean position of species from the global pool (present in all assemblages).

  • FNND Functional Mean Nearest Neighbour Distance: the weighted distance to the nearest neighbor within the assemblage

  • FMPD Functional Mean Pairwise Distance: the mean weighted distance between all species pairs.

Code
# Data Preparation ----
station_morpho <- morpho_data %>% 
  select(station) %>% 
  filter(!station %in% c("F0315", "F0367", "X0568")) %>%  
  arrange(station) %>%
  distinct()

station_sp <- rbind(data_biomass_2002_2019, data_biomass_2021, data_biomass_2022) %>%
  as.data.frame() %>%
  rename("species" = "Nom_Scientifique",
         "station" = "Code_Station") %>%
  left_join(metadata) %>%
  select(species, Tot_V_HV, volume_filtered, station) %>%
  # Divide biomass by the volume filtered at each trawl (g.m3)
  mutate(biomass_cpu = (Tot_V_HV / volume_filtered) * 1000) %>%
  select(species, biomass_cpu, station) %>%
  replace(is.na(.), 0) %>%
  group_by(species, station) %>%
  mutate(biomass = sum(biomass_cpu)) %>%
  select(-biomass_cpu) %>%
  distinct() %>%
  tidyr::pivot_wider(names_from = species, values_from = biomass) %>%
  replace(is.na(.), 0) %>%
  arrange(station) %>%
  filter(!station%in%c("H0411","L0731", "L0736")) %>%
  tibble::column_to_rownames(var = "station") %>%
  # Change species names
  rename("Lampanyctus_ater" = "Nannobrachium_atrum") %>%
  rename("Cyclothone_sp" = "Cyclothone") %>%
  rename("Stomias_boa" = "Stomias_boa_boa") %>%
  select(order(colnames(.))) %>%
  as.matrix()

# if we want only presence-absence 
#station_sp <- replace(station_sp, station_sp > 0, 1)

# dist matrix (gower)
sp_dist_fish <- mFD::funct.dist(
  sp_tr         = fish_traits,
  tr_cat        = fish_traits_cat,
  metric        = "gower",
  scale_euclid  = "scale_center",
  ordinal_var   = "classic",
  weight_type   = "equal",
  stop_if_NA    = TRUE)

# Compute functional spaces 
fspaces_quality_fish <- mFD::quality.fspaces(
  sp_dist             = sp_dist_fish,
  maxdim_pcoa         = 10,
  deviation_weighting = "absolute",
  fdist_scaling       = FALSE,
  fdendro             = "average")

# PC coord
sp_faxes_coord_fish <- fspaces_quality_fish$"details_fspaces"$"sp_pc_coord"

# Calculate functional diversity for the observed data ----
obsFD <- mFD::alpha.fd.multidim(
  sp_faxes_coord = sp_faxes_coord_fish[, c("PC1", "PC2", "PC3", "PC4")],
  asb_sp_w = station_sp,
  scaling = TRUE,
  check_input = TRUE,
  details_returned = F
)

obsFD_div <- obsFD$functional_diversity_indices

# Null model ----
# Define the number of replications
nb_rep <- 10

# Initialize a list to store results of random functional diversity calculations for each index
indices_names <- colnames(obsFD_div)
resultsRandomFD <- list()

for (index_name in indices_names) {
  resultsRandomFD[[index_name]] <- matrix(
    NA,
    nrow = nrow(station_sp),
    ncol = nb_rep,
    dimnames = list(rownames(station_sp), paste0("Sim.", 1:nb_rep))
  )
}

# Perform randomization and calculate functional diversity for each replication
for (rep in 1:nb_rep) {
  randomize_mx <- picante::randomizeMatrix(
    samp = station_sp,
    null.model = "frequency", # richness, independentswap, trialswap
    iterations = 10
  )
  
  simFD_cal <- mFD::alpha.fd.multidim(
    sp_faxes_coord = sp_faxes_coord_fish[, c("PC1", "PC2", "PC3", "PC4")],
    asb_sp_w = randomize_mx,
    scaling = TRUE,
    check_input = TRUE,
    details_returned = F
  )
  
  simFD_div <- simFD_cal$functional_diversity_indices
  
  for (index_name in indices_names) {
    simFD_index <- simFD_div[, index_name]
    
    # Ensure that simFD_index has the same length as the number of rows in station_sp
    if(length(simFD_index) == nrow(station_sp)) {
      resultsRandomFD[[index_name]][, rep] <- simFD_index
    } else {
      stop(paste("The length of", index_name, "does not match the number of rows in station_sp"))
    }
  }
}

# Initialize dataframes to store mean, standard deviation, effect size, and standardized effect size
meanNullFD <- data.frame(matrix(NA, nrow = nrow(station_sp), ncol = length(indices_names)))
sdNullFD <- data.frame(matrix(NA, nrow = nrow(station_sp), ncol = length(indices_names)))
ES_FD <- data.frame(matrix(NA, nrow = nrow(station_sp), ncol = length(indices_names)))
SES_FD <- data.frame(matrix(NA, nrow = nrow(station_sp), ncol = length(indices_names)))

# Set column names for the dataframes
colnames(meanNullFD) <- indices_names
colnames(sdNullFD) <- indices_names
colnames(ES_FD) <- indices_names
colnames(SES_FD) <- indices_names

# Calculate statistics for each index
for (index_name in indices_names) {
  # Calculate mean and standard deviation of null model FD values for each index
  meanNullFD[, index_name] <- rowMeans(resultsRandomFD[[index_name]], na.rm = TRUE)
  sdNullFD[, index_name] <- apply(resultsRandomFD[[index_name]], 1, sd, na.rm = TRUE)
  
  # Calculate effect size and standardized effect size for each index
  ES_FD[, index_name] <- obsFD_div[, index_name] - meanNullFD[, index_name]
  SES_FD[, index_name] <- ES_FD[, index_name] / sdNullFD[, index_name]
}

# Combine all results into a single dataframe
results_df <- cbind(
  # obsFD_div,
  # meanNullFD = meanNullFD,
  # sdNullFD = sdNullFD,
  # ES_FD = ES_FD,
  SES_FD = SES_FD
)

# Add row names
rownames(results_df) <- rownames(station_sp)

# Plot  ----
# Output the results dataframe
results_df_plot <- results_df %>%
  tibble::rownames_to_column(var = "station") %>% 
  inner_join(metadata %>% select(station, depth), by = "station") %>%
  mutate(
    depth_layer = case_when(
      between(depth, 0, 174) ~ "Epipelagic",
      between(depth, 175, 699) ~ "Upper mesopelagic",
      between(depth, 700, 999) ~ "Lower mesopelagic",
      between(depth, 1000, 2000) ~ "Bathypelagic"
    )
  ) %>% 
  select(-station) %>%
  tidyr::pivot_longer(!c(depth, depth_layer),
                      values_to = "values",
                      names_to = "indice") %>%
  mutate(indice = stringr::str_replace(indice, "^SES_", "")) %>%
  filter(
    indice %in% c(
      "FD.fric",
      "FD.fdis",
      "FD.fmpd",
      "FD.fnnd",
      "FD.feve",
      "FD.fdiv",
      "FD.fori",
      "FD.fspe"
    )
  )

results_df_plot$depth_layer <- factor(results_df_plot$depth_layer, 
                                      levels = c("Epipelagic", "Upper mesopelagic",
                                                 "Lower mesopelagic", "Bathypelagic"))
results_df_plot$indice <- factor(results_df_plot$indice,
                                 levels = c("FD.fric",
                                            "FD.fdis",
                                            "FD.fdiv",
                                            "FD.feve",
                                            "FD.fori",
                                            "FD.fspe",
                                            "FD.fnnd",
                                            "FD.fmpd"
                                            ))

ggplot(results_df_plot, aes(x = depth_layer, y = values, fill = depth_layer)) +
  facet_wrap(~indice, scales = "free") +
  geom_point(alpha = 0.5, size = 1, position = position_jitter(width = 0.2), aes(col=depth_layer)) +
  geom_boxplot(alpha = 0.1, outlier.shape = NA, width = 0.5, aes(col=depth_layer, fill=depth_layer))+
  scale_color_manual(values = c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +
  scale_fill_manual(values = c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +
  labs(
    x = "",
    y = "Standard Effect Size (SES)"
  ) +
  theme_light() +
  theme(axis.text.x = element_blank(),
        strip.text.x = element_text(
          size = 13,
          face = "bold",
          color = "gray20"
        ),
        strip.background = element_rect(fill = "white"),
        axis.title = element_text(size = 13),
        axis.text = element_text(size = 13))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "black", linewidth = 0.8)+
  guides(col="none", fill="none")

Isotopes data

  • 4 sp manquantes en d15N = P. coregonoides, E. pelacanoides, H. anomala, S. bathyphilus

  • functional space : colours according to d15N values

  • significant relationships between d15n and traits

  • categorical traits

Plan

I. Key traits and foraging strategies by depth layer

Figure 1: Lateral-view photographs

Figure 2: Functional spaces

  • At the community level & by depth layer

  • Only PC1-PC2 and PC3-PC4 in appendices?

Table 1: significant traits

  • which traits significantly influence each axis?

  • at the community level & by depth layer

Figure 3: Community-Weighted Mean

II. Functional diversity indices

Figure 4: Standardized Effect Size

  • Trait convergence in the epipelagic layer and trait divergence in the bathypelagic layer

III. Functional rarity

Figure 5: Uniqueness-geographical restrectiveness